###SETUP####
rm(list = ls())
library(tidyverse) 
library(plyr)
library(dplyr)
library(cregg)
library(ggplot2)
library(gridExtra)
library(grid)
library(stringr)

# to replicate, set the working directory to correct path on your computer
setwd("your.path.here")

dat <- read_csv(file ="HKM_profiledata.csv")
dat2 <- read_csv(file ="HKM_respondentdata.csv")

####CLEANING & VARIABLE CONSTRUCTION####

#omitting missing values
fdata <- na.omit(dat)

#Recode Gender Ideology Scale & Fairness Scale to Numeric
fdata <-
  fdata %>%
  mutate(
modernsexism_1n = recode(modernsexism_1,
                         "Strongly Disagree" = 0, "Somewhat Disagree" = 1,
                         "Neither agree nor disagree" = 2, "Somewhat Agree" = 3,
                         "Strongly Agree"= 4),
modernsexism_2n = recode(modernsexism_2,
                         "Strongly Disagree" = 0, "Somewhat Disagree" = 1,
                         "Neither agree nor disagree" = 2, "Somewhat Agree" = 3,
                         "Strongly Agree"= 4),
modernsexism_3n = recode(modernsexism_3,
                         "Strongly Disagree" = 4, "Somewhat Disagree" = 3,
                         "Neither agree nor disagree" = 2, "Somewhat Agree" = 1,
                         "Strongly Agree"= 0),
modernsexism_4n = recode(modernsexism_4,
                         "Strongly Disagree" = 0, "Somewhat Disagree" = 1,
                         "Neither agree nor disagree" = 2, "Somewhat Agree" = 3,
                         "Strongly Agree"= 4),
profile_fairnum = recode(profile_fair,
                         "Very unfair" = 0, "Somewhat unfair" = 1,
                         "Neither fair nor unfair" = 2, 
                         "Somewhat fair" = 3,
                         "Very fair" = 4)
          )

#Standardize Fairness Outcome Measure

fdata <- fdata %>% mutate_each_(funs(scale(.) %>% as.vector), 
                                vars=c("profile_fairnum"))

#Generate Modern Sexism Measure
fdata$modernsexism_sum=rowSums(cbind(fdata$modernsexism_1n,fdata$modernsexism_2n,fdata$modernsexism_3n,fdata$modernsexism_4n),na.rm=TRUE)
fdata$sexismscale[fdata$modernsexism_sum <= median(fdata$modernsexism_sum)] <- 1L
fdata$sexismscale[fdata$modernsexism_sum > median(fdata$modernsexism_sum)] <- 2L
fdata$sexismscale<- factor(fdata$sexismscale, 1:2, c("Low Sexism", "High Sexism"))

#Convert attribute variables to factors
fdata$childcare=as.factor(fdata$childcare_level)
fdata$contribution=as.factor(fdata$contribution_level)
fdata$fearnings=as.factor(fdata$fearnings_level)
fdata$mearnings=as.factor(fdata$mearnings_level)
fdata$division=as.factor(fdata$division_level)
fdata$dem_gender=as.factor(fdata$dem_gender)
fdata$sexismscale_factor=as.factor(fdata$sexismscale)

#Numeric gender variable for balance analysis
fdata$female=ifelse(fdata$dem_gender=='Female', 1,0)

#Set level labels
levels(fdata$childcare) <- c("At least 1", "None")
levels(fdata$fearnings) <- c ("Wife Unemployed", "Wife Works, earns 100k", "Wife Works, earns 50k")
levels(fdata$mearnings) <- c ("Husband Unemployed", "Husband Works, earns 100k", "Husband Works, earns 50k")
levels(fdata$division) <- c ("Wife spends about same time", "Wife spends much less time", "Wife spends much more time")
levels(fdata$contribution) <- c ("Husband rarely contributes", "Husband regularly contributes", "Husband sometimes contributes")
levels(fdata$sexismscale_factor) <- c ("Low", "High")

#Reorder factors
fdata$childcare <- factor(fdata$childcare,levels(fdata$childcare)[c(2,1)])
fdata$fearnings <- factor(fdata$fearnings,levels(fdata$fearnings)[c(1,3,2)])
fdata$mearnings <- factor(fdata$mearnings,levels(fdata$mearnings)[c(1,3,2)])
fdata$contribution <- factor(fdata$contribution,levels(fdata$contribution)[c(1,3,2)])

#Set Baselines
baselines <- list()
baselines$childcare <- "None"
baselines$fearnings <- "Wife Unemployed"
baselines$mearnings <- "Husband Unemployed"
baselines$division <- "Wife spends about same time"
baselines$contribution <-"Husand rarely contributes"

####DEFINE FORMULAS####

f1<- profile1_chosen ~ childcare + fearnings + mearnings + 
  division + contribution

f2<- profile_fairnum ~ childcare + fearnings + mearnings + 
  division + contribution

f3<-profile1_chosen ~ fearnings + mearnings + 
  division + contribution

f4<-profile_fairnum ~ fearnings + mearnings + 
  division + contribution

f5<- female ~ childcare + fearnings + mearnings + 
  division + contribution

f6<- modernsexism_sum ~ childcare + fearnings + mearnings + 
  division + contribution

####ESTIMATE MODELS####

#AMCE (Outcome=Profile Choice)
amce1 <- cj(fdata, 
           f1,
            id = ~rid, estimate = "amce")

#AMCE (Outcome=Fairness Perception)
amce2 <- cj(fdata, 
            f2,
            id = ~rid, estimate = "amce")

#Conditional AMCE by Gender (Outcome=Profile Choice)
amces_gender <- cj(fdata, 
                   f1,
                   id = ~rid, estimate = "amce",
                   by = ~dem_gender)

amces_genderdiff <- cj(fdata, 
                       f1,
                       id = ~rid, estimate = "amce_diff",
                       by = ~dem_gender)

#Conditional AMCE by Gender (Outcome=Fairness Perception)
amces_gender2 <- cj(fdata, 
                    f2,
                    id = ~rid, estimate = "amce",
                    by = ~dem_gender)

amces_genderdiff2 <- cj(fdata, 
                        f2,
                        id = ~rid, estimate = "amce_diff",
                        by = ~dem_gender)

#Conditional AMCE by Gender Ideology (Outcome=Profile Choice)
amces_sexism <- cj(fdata,
                   f1,
                   id = ~rid, estimate = "amce",
                   by = ~sexismscale_factor)

amces_sexismdiff <- cj(fdata,
                       f1,
                       id = ~rid, estimate = "amce_diff",
                       by = ~sexismscale_factor)

#Conditional AMCE by Gender Ideology (Outcome=Fairness Perception)
amces_sexism2 <- cj(fdata, 
                    f2,
                    id = ~rid, estimate = "amce",
                    by = ~sexismscale_factor)

amces_sexismdiff2 <- cj(fdata, 
                        f2,
                        id = ~rid, estimate = "amce_diff",
                        by = ~sexismscale_factor)

#ACIE by Childcare Availability (Outcome=Profile Choice)
amces_childcare <- cj(fdata, 
                      f3,
                      id = ~rid, estimate = "amce",
                      by = ~childcare)

amces_childcarediff <- cj(fdata, 
                          f3,
                          id = ~rid, estimate = "amce_diff",
                          by = ~childcare)

#ACIE by Childcare Availability (Outcome=Fairness Perception)
amces_childcare2 <- cj(fdata, 
                       f4,
                       id = ~rid, estimate = "amce",
                       by = ~childcare)

amces_childcarediff2 <- cj(fdata, 
                           f4,
                           id = ~rid, estimate = "amce_diff",
                           by = ~childcare)

#DIAGNOSTICS

#Frequency of attribute values

att_freq<-cj_freqs(fdata, f1, id = ~rid)

#Attribute Balance across covariates (Gender)

att_balance1<-cj(fdata,
                 f5,
                 id = ~rid, estimate = "mm",
                )

#Attribute Balance across covariates (Gender Ideology)

att_balance2<-cj(fdata,
                 f6,
                 id = ~rid, estimate = "mm",
)

####FORMAT & OUTPUT####

#FIGURE 1 
p1<-plot(amce1) + ggplot2::theme(legend.position = "none") + scale_color_manual(values = rep("black", 5)) + ggtitle("Profile Choice") +
  scale_y_discrete(labels = c("(childcare)"= "CHILDCARE OPTIONS", 
                              "(fearnings)" = "WIFE WORK STATUS", 
                              "Wife Unemployed" = "Unemployed",
                              "Wife Works, earns 100k" = "Works, earns 100k", 
                              "Wife Works, earns 50k" = "Works, earns 50k", 
                              "(mearnings)" = "HUSBAND WORK STATUS",
                              "Husband Unemployed" = "Unemployed",
                              "Husband Works, earns 100k" = "Works, earns 100k", 
                              "Husband Works, earns 50k" = "Works, earns 50k",
                              "(division)" = "TIME DIVISION ON HH CHORES",
                              "(contribution)" = "CONTRIBUTION TO FEMININE TASKS"))


p2<-plot(amce2) + ggplot2::theme(legend.position = "none",
                                 axis.text.y = element_blank()) + scale_color_manual(values = rep("black", 5)) + ggtitle("Fairness Evaluation")

g <- arrangeGrob(p1, p2, widths = c(1.7,1), nrow=1)
ggsave(file="tabsfigs/fig1.png", g)

#FIGURE 2

p3<-plot(rbind(amces_gender, amces_genderdiff)) + ggplot2::facet_wrap(~BY, ncol = 3L) + 
  theme(legend.position = "none") + scale_color_manual(values = rep("black", 5)) + 
  ggtitle("Profile Choice, by Gender") +
  scale_y_discrete(labels = c("(childcare)"= "CHILDCARE OPTIONS", 
                              "(fearnings)" = "WIFE WORK STATUS", 
                              "Wife Unemployed" = "Unemployed",
                              "Wife Works, earns 100k" = "Works, earns 100k", 
                              "Wife Works, earns 50k" = "Works, earns 50k", 
                              "(mearnings)" = "HUSBAND WORK STATUS",
                              "Husband Unemployed" = "Unemployed",
                              "Husband Works, earns 100k" = "Works, earns 100k", 
                              "Husband Works, earns 50k" = "Works, earns 50k",
                              "(division)" = "TIME DIVISION ON HH CHORES",
                              "(contribution)" = "CONTRIBUTION TO FEMININE TASKS"))

ggsave(file="tabsfigs/fig2.png", p3)


#FIGURE 3

p4<-plot(rbind(amces_gender2, amces_genderdiff2)) + ggplot2::facet_wrap(~BY, ncol = 3L) + 
  theme(legend.position = "none") + scale_color_manual(values = rep("black", 5)) + 
  ggtitle("Fairness Evaluation, by Gender") +
  scale_y_discrete(labels = c("(childcare)"= "CHILDCARE OPTIONS", 
                              "(fearnings)" = "WIFE WORK STATUS", 
                              "Wife Unemployed" = "Unemployed",
                              "Wife Works, earns 100k" = "Works, earns 100k", 
                              "Wife Works, earns 50k" = "Works, earns 50k", 
                              "(mearnings)" = "HUSBAND WORK STATUS",
                              "Husband Unemployed" = "Unemployed",
                              "Husband Works, earns 100k" = "Works, earns 100k", 
                              "Husband Works, earns 50k" = "Works, earns 50k",
                              "(division)" = "TIME DIVISION ON HH CHORES",
                              "(contribution)" = "CONTRIBUTION TO FEMININE TASKS"))

ggsave(file="tabsfigs/fig3.png", p4)


#FIGURE 4

p5<-plot(rbind(amces_sexism, amces_sexismdiff)) + ggplot2::facet_wrap(~BY, ncol = 3L) + 
  theme(legend.position = "none") + scale_color_manual(values = rep("black", 5)) + 
  ggtitle("Profile Choice, by Sexism Level") +
  scale_y_discrete(labels = c("(childcare)"= "CHILDCARE OPTIONS", 
                              "(fearnings)" = "WIFE WORK STATUS", 
                              "Wife Unemployed" = "Unemployed",
                              "Wife Works, earns 100k" = "Works, earns 100k", 
                              "Wife Works, earns 50k" = "Works, earns 50k", 
                              "(mearnings)" = "HUSBAND WORK STATUS",
                              "Husband Unemployed" = "Unemployed",
                              "Husband Works, earns 100k" = "Works, earns 100k", 
                              "Husband Works, earns 50k" = "Works, earns 50k",
                              "(division)" = "TIME DIVISION ON HH CHORES",
                              "(contribution)" = "CONTRIBUTION TO FEMININE TASKS"))


ggsave(file="tabsfigs/fig4.png", p5)

# FIGURE 5

p6<-plot(rbind(amces_sexism2, amces_sexismdiff2)) + ggplot2::facet_wrap(~BY, ncol = 3L) + 
  theme(legend.position = "none") + scale_color_manual(values = rep("black", 5)) + 
  ggtitle("Fairness Evaluation, by Sexism Level") +
  scale_y_discrete(labels = c("(childcare)"= "CHILDCARE OPTIONS", 
                              "(fearnings)" = "WIFE WORK STATUS", 
                              "Wife Unemployed" = "Unemployed",
                              "Wife Works, earns 100k" = "Works, earns 100k", 
                              "Wife Works, earns 50k" = "Works, earns 50k", 
                              "(mearnings)" = "HUSBAND WORK STATUS",
                              "Husband Unemployed" = "Unemployed",
                              "Husband Works, earns 100k" = "Works, earns 100k", 
                              "Husband Works, earns 50k" = "Works, earns 50k",
                              "(division)" = "TIME DIVISION ON HH CHORES",
                              "(contribution)" = "CONTRIBUTION TO FEMININE TASKS"))


ggsave(file="tabsfigs/fig5.png", p6)

# FIGURE 6

p7<-plot(rbind(amces_childcare, amces_childcarediff)) + ggplot2::facet_wrap(~BY, ncol = 3L) + 
  theme(legend.position = "none") + scale_color_manual(values = rep("black", 5)) + 
  ggtitle("Profile Choice, by Childcare Availability") +
  scale_y_discrete(labels = c( "(fearnings)" = "WIFE WORK STATUS", 
                              "Wife Unemployed" = "Unemployed",
                              "Wife Works, earns 100k" = "Works, earns 100k", 
                              "Wife Works, earns 50k" = "Works, earns 50k", 
                              "(mearnings)" = "HUSBAND WORK STATUS",
                              "Husband Unemployed" = "Unemployed",
                              "Husband Works, earns 100k" = "Works, earns 100k", 
                              "Husband Works, earns 50k" = "Works, earns 50k",
                              "(division)" = "TIME DIVISION ON HH CHORES",
                              "(contribution)" = "CONTRIBUTION TO FEMININE TASKS"))

ggsave(file="tabsfigs/fig6.png", p7)

# FIGURE 7

p8<-plot(rbind(amces_childcare2, amces_childcarediff2)) + ggplot2::facet_wrap(~BY, ncol = 3L) + 
  theme(legend.position = "none") + scale_color_manual(values = rep("black", 5)) + 
  ggtitle("Fairness Evaluation, by Childcare Availability") +
  scale_y_discrete(labels = c( "(fearnings)" = "WIFE WORK STATUS", 
                               "Wife Unemployed" = "Unemployed",
                               "Wife Works, earns 100k" = "Works, earns 100k", 
                               "Wife Works, earns 50k" = "Works, earns 50k", 
                               "(mearnings)" = "HUSBAND WORK STATUS",
                               "Husband Unemployed" = "Unemployed",
                               "Husband Works, earns 100k" = "Works, earns 100k", 
                               "Husband Works, earns 50k" = "Works, earns 50k",
                               "(division)" = "TIME DIVISION ON HH CHORES",
                               "(contribution)" = "CONTRIBUTION TO FEMININE TASKS"))


ggsave(file="tabsfigs/fig7.png", p8)


## FIGURE A1

p9<- plot(att_freq) + ggplot2 :: theme(legend.position="none") + scale_fill_grey() + ggtitle("Frequency of Attributes") + scale_x_discrete(labels = c("(childcare)"= "CHILDCARE OPTIONS", 
                              "(fearnings)" = "WIFE WORK STATUS", 
                              "Wife Unemployed" = "Unemployed",
                              "Wife Works, earns 100k" = "Works, earns 100k", 
                              "Wife Works, earns 50k" = "Works, earns 50k", 
                              "(mearnings)" = "HUSBAND WORK STATUS",
                              "Husband Unemployed" = "Unemployed",
                              "Husband Works, earns 100k" = "Works, earns 100k", 
                              "Husband Works, earns 50k" = "Works, earns 50k",
                              "(division)" = "TIME DIVISION ON HH CHORES",
                              "(contribution)" = "CONTRIBUTION TO FEMININE TASKS"))


ggsave(file="tabsfigs/figA1.png", p9)

# FIGURE A2

p10 <- plot(att_balance1, vline=mean(fdata$female)) + ggplot2:: theme(legend.position = "none") + scale_color_manual(values = rep("black", 5)) + 
  ggtitle("Gender (Female=1)") +
  scale_y_discrete(labels = c("(childcare)"= "CHILDCARE OPTIONS", 
                              "(fearnings)" = "WIFE WORK STATUS", 
                              "Wife Unemployed" = "Unemployed",
                              "Wife Works, earns 100k" = "Works, earns 100k", 
                              "Wife Works, earns 50k" = "Works, earns 50k", 
                              "(mearnings)" = "HUSBAND WORK STATUS",
                              "Husband Unemployed" = "Unemployed",
                              "Husband Works, earns 100k" = "Works, earns 100k", 
                              "Husband Works, earns 50k" = "Works, earns 50k",
                              "(division)" = "TIME DIVISION ON HH CHORES",
                              "(contribution)" = "CONTRIBUTION TO FEMININE TASKS"))
                             

ggsave(file="tabsfigs/figA2.png", p10)

## FIGURE A3

p11 <- plot(att_balance2, vline=mean(fdata$modernsexism_sum)) + ggplot2:: theme(legend.position = "none") + scale_color_manual(values = rep("black", 5)) + 
  ggtitle("Gender Ideology (Sexism Scale)") +
  scale_y_discrete(labels = c("(childcare)"= "CHILDCARE OPTIONS", 
                              "(fearnings)" = "WIFE WORK STATUS", 
                              "Wife Unemployed" = "Unemployed",
                              "Wife Works, earns 100k" = "Works, earns 100k", 
                              "Wife Works, earns 50k" = "Works, earns 50k", 
                              "(mearnings)" = "HUSBAND WORK STATUS",
                              "Husband Unemployed" = "Unemployed",
                              "Husband Works, earns 100k" = "Works, earns 100k", 
                              "Husband Works, earns 50k" = "Works, earns 50k",
                              "(division)" = "TIME DIVISION ON HH CHORES",
                              "(contribution)" = "CONTRIBUTION TO FEMININE TASKS"))


ggsave(file="tabsfigs/figA3.png", p11)

## TABLE A1: Respondent Summary Stats
# Created manually from output


#GENDER
prop.table(table(dat2$dem_gender))

#AGE
mean(dat2$age_consent)

#EDUCATION
dat2$bachelors<-ifelse((dat2$edu=="Bachelors degree"|
                        dat2$edu=="Masters degree"|
                        dat2$edu=="Professional degree (ex: JD, MBA, MD)"|
                        dat2$edu=="Doctorate"),1,0)
mean(dat2$bachelors)

#RACE & ETHNICITY
dat2$white<-str_detect(dat2$race, "White")
dat2$black<-str_detect(dat2$race, "Black or African American")
dat2$hispanic<-str_detect(dat2$race, "Hispanic or Latino")
dat2$asian<- str_detect(dat2$race, "Asian or Asian American")
dat2$native<- str_detect(dat2$race, "Native American")
dat2$other<- str_detect(dat2$race, "Other")

prop.table(table(dat2$white))
prop.table(table(dat2$black))
prop.table(table(dat2$hispanic))
prop.table(table(dat2$asian))
prop.table(table(dat2$native))
prop.table(table(dat2$other))

#HH INCOME
prop.table(table(dat2$income_house))

dat2$income1<-ifelse((dat2$income_house=="Under $10,000"|
                     dat2$income_house=="$10,000 to $20,000"),1,0)

dat2$income2<-ifelse((dat2$income_house=="$20,000 to $30,000"|
                     dat2$income_house=="$30,000 to $40,000"|
                     dat2$income_house=="$40,000 to $50,000" ),1,0)

dat2$income3<-ifelse((dat2$income_house=="$50,000 to $60,000"|
                        dat2$income_house=="$60,000 to $70,000"|
                        dat2$income_house=="$70,000 to $80,000" |
                        dat2$income_house=="$80,000 to $90,000" |
                        dat2$income_house=="$90,000 to $100,000"),1,0)

dat2$income4<-ifelse((dat2$income_house=="$100,000 to $150,000"),1,0)

dat2$income5<-ifelse((dat2$income_house=="Over $150,000"),1,0)

prop.table(table(dat2$income1))
prop.table(table(dat2$income2))
prop.table(table(dat2$income3))
prop.table(table(dat2$income4))
prop.table(table(dat2$income5))

##EMPLOYMENT
prop.table(table(dat2$employ))

##PARTY ID
prop.table(table(dat2$pid))
